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Abstract 

We propose a first-principles time-dependent density functional theoretical (TDDFT) approach in mo- 
mentum (V) space for quantitative study of electron transport in molecular devices under arbitrary biases. 
In this approach, the basic equation of motion is a time-dependent integrodifferential equation obtained by 
Fourier transform of the time-dependent Kohn-Sham (TDKS) equation in spatial coordinate (TZ) space. It 
is formally exact and includes all the effects and information of the electron transport in molecular devices. 
The electron wavefunction is calculated by solving this equation in a closed finite P-space volume. This 
approach is free of self-energy function and memory term and beyond the wide-band limit ( WBL) . The fea- 
sibility and power of the approach are demonstrated by the calculation of current through one-dimensional 
(ID) systems. 

PACS numbers: 73.63.-b,85.65.-Hh,71.15.Pd,31. 15.ee 
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Building molecular electronic devices using individual molecules is one of the ultimate goals in 
nanotechnology . A typical molecular device consists of a molecule coupled to two (or 

more) electrodes [7|] . Electrons in the molecular device can go anywhere along the electrodes during 



the electron transport and electron wavefunctions may extend to infinity in spatial coordinate (7^) 
space. Thus, accurate calculation of the electron wavefunctions in the TZ space is an intractable 
task for the study of electron transport. 

To resolve this problem one of the widely used methods is to separate the molecule from the 



electrodes and treat it as an open system. More specifical 



left (L), central (C), and right (R) zones 



S 



Id, 



111, 



y, t he molecular device is partitioned into 



3,Q,Q. 



The C zone is chosen to include 



the molecule and some atomic layers of the electrodes so that the Hamiltonian and electron density 
of the L and R zones are accurately described by the equilibrium bulk ones before the bias applied 
Q]. To calculate the electron wavefunctions, one separates the C zone from the L and R zones and 
treat the C zone as an open system. The effect of the L and R zon es ( environmeiit) on the C zone 



0,y,Q,Q,[3,Q, 



14 



15|, 



la ] . After taking 



(open system) is characterized by a self-energy function 
this effect into account, an equation of motion for the C zone is derived. Solving this equation 
one achieves the electron wavefunctions for further calculation. This scheme has been extensively 



0,y,Q, 
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i7i, Q,y. 



applied to the calculation of steady-state current in molecular devices 

However, this scheme is not completely applicable when a time-dependent bias is applied. In 
this case, the equation of motion of the C zone is derived from the time-dependent Kohn-Sham 



(TDKS) equationJlill3jl4[. 



a memory term 



14 



151, 



t not only contains a two-time self-energy function but also includes 



16l |. The self-energy function is related to a two-time Green function 



that can only be calculated under the wide-band limit (WBL) [20], an approximation that neglects 
the energy dependence of the coupling between the molecule and electrodes [l^. Beyond the WBL, 
the two-time Green function is governed by a double integral equation which is generally unsolvable 
[isl . The memory term depends on the two-time self-energy function. To calculate the memory 
term at time t, one has to recalculate and store the self-energy function at all the past time steps 
before this time. As a result, computer resources such as CPU time and random-access memory 
(RAM) required in the calculation become increasingly and thus extremely large with increase of 
number of time steps. 



To overcome these difficu^ 
solving TDKS equation 
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ties, a computationally feasible scheme was proposed based on directly 
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IJ]. In this scheme, the TDKS equation is discretized in the whole 
TZ space and converted to a matrix equation first, then an equation of motion is derived from 
the matrix equation for the C zone, and finally the equation of motion is solved to calculate the 



2 



electron wavefunctions. This scheme has been successfully applied to the study of ID model systems 
12,ll3, lljl. The time-dependent currents achieved tend to the steady-state currents obtained from 

nn n 

the Landauer formula [Ij, ll3|] and Floquet method [Ij] after a long time, providing a benchmark 
for the currents through the ID model systems. The most striking characteristic of this pioneering 
scheme is that the effect of the L and R zones on the C-zone electrons and the reflection of the 
time-dependent wavefunction on the boundaries are completely taken good care of by a transparent 
boundary condition. This scheme can indeed eliminate the explicit dependence of the equation of 
motion on the two-time self-energy function. However, the equation of motion still contains a 
memory term used to represent the effect of the L and R zones and remove the reflection of the 
electron wavefunction on the boundaries. As a consequence, the computer resources required in 
the calculation also increase with the number of time steps. 

In fact, the difficulties encountered in the 7^-space calculations at least come from one of the 
memory term and self-energy function used to characterize the effect of the L and R zones on the 
electrons in the C zone. Since the electron wavefunction can extend to the infinity in the molecular 
device during the transport, the memory term or/and self-energy function can not be disregarded 
no matter how large the C zone is chosen. Because of this reason, quantitative study of electron 
transient transport in the TZ space is currently a computational challenge for realistic molecular 
devices. 

However, the momentum of electron is always finite and less than certain maximum value 
kjnax in the molecular device. The probability of electron with momentum greater than k^ax is 
negligible or zero. Thus, in momentum (V) space the electron wavefunction is localized (e.g., the 
wavefunction of a free electron spreads out to the whole space in the TZ space but localizes in a very 
small volume depicted by a 5 function in the V space) [21] and can be calculated in a finite volume 
as long as the boundary of the volume is set at the place with properly large momentum kj^ax* 

In 

this case, all the troublesome terms related to the open system in the TZ space, such as the self- 
energy function and memory term, will be wiped out completely and the difficulties encountered 
in the 7^-space calculations can be resolved in the V space. Based on this idea, we propose in this 
work a first-principles time-dependent density functional theoretical (TDDFT) approach in the V 
space for the quantitative study of electron transient transport in realistic molecular devices under 
arbitrary biases. 

The molecular device is in an equilibrium state described by a unique temperature and chemical 
potential when time t < 0. The charge of the two electrodes is perfectly balanced and no current 
flows through the device. At time t = a bias of voltage Vb (t) is applied to the electrodes. This 
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bias drives the molecular device out of equilibrium and induces the current through the device. 

In the TZ space, the time-dependent electron wavefunction governed by the TDKS 

equation within TDDFT (atomic units are used throughout the paper) 

i^V(r, t) = [ho (r) + VD (r, t)] i,{r, t), (1) 

where, ho (r) is the unperturbed KS Hamiltonian of the electron and (r, t) is the driving potential 
characterizing the effect of bias on the electron. Within the partition above, the unperturbed KS 
Hamiltonian can be written as /iq (r) = — + vq (r), where, vq (r) = vbko (r) with a = L 
and R for the L and R zones and vq (r) = Vext (r) + ^e//(r,0) for the C zone. Here, vbko (r) 
is the bulk potential of the electrodes, Vext{^) is the external potential related to the interaction 



between the electron and nuclei of the system 



22], and fe//(r, 0) is the effective potential of the 



electron comprising the Hartree potential vh (r, 0) and exchange-correlation potential fxc(r, 0). For 
the metallic electrodes considered here, the driving potential vd (r, t) is represented by (t) and 
vji (t) in the L and R zones, respectively, and vj^ (t) — vr (t) = Vh (t). These potentials undergo 
uniform time-dependent shifts due to the screening effect and thus do not chan ge w ith spatial 
coordinates if Vh{t) is slowly varying during a typical time scale (less than a fs) 1^. In the C 
zone, the driving potential vd (r, t) is represented by vc (r, t) and vc (r, t) = Veff{r, t) — Veff{r, 0) 
is the perturbing potential induced by the field of bias, where Veff{r, t) is the bias-induced effective 
potential including the Hartree potential vh i^,t) and exchange-correlation potential Vxci^,t) 23l |. 
The bias-induced effective potential is associated with the response of the C-zone electrons to the 
field of bias. It may change dramatically and nonlinearly with spatial coordinates and can only be 
calculated self-consistently with the C-zone electron wavef unctions. Theprocedure proposed here 



231]. 



is nonperturbative and goes beyond the linear response approximation 

The relation between the time-dependent 7?.-space electron wavefunction ip{r,t) and 
"P-space electron wavefunction (j)(k,t) is given by the Fourier transform iplrji) = 
(27r)~^/^ J dk(/)(k, t) exp (zk • r). Applying this transform to Eq. ([T]), we obtain the V-space TDKS 
integrodifferential equation 



i-(/)(k,t) = J dk' [Ho{k,k') + VDiKk',t)] Hk',t), (2) 

where, Hq (k, k') = (k — k') /2 + Vq (k, k') is the T'-space unperturbed Hamiltonian, Vq (k, k') 
and Voik, k', t) are the P-space potentials calculated from the 7^-space potentials vq (r) and vd (r, t) 
by the transform V^(k, k') = (27r)^^ / dvva{v) exp [i (k! — k) • r] with a = and D, respectively. 
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To calculate the 7^-space electron wavefunction, we solve Eq. ([2|) numerically in a finite 3D 
volume [—ka max — ^Q. — ^ «max|Q^ — x,y,z\. The integrals in Eg . (121) a re discretized using 
the generalized Legendre-Gauss-Lobatto pseudospectral method [22, \2S, l25| on 3D grid points 



[0, 1, • • • ,Na\a = x,y,z] in the finite 3D volume. Since the 7^-space electron wavefunction is 
localized the boundary condition for solving Eq. ^ is simply that the electron wavefunction 
is zero on the boundary. After discretization, Eq. ^ is converted to a time-dependent matrix 
equation 

i^*(t) = [Ho+Vz)(t)]$(t), (3) 

where, is an A'^c-dimensional vector with the components being the P-space electron wave- 
function (/>(k, t) on the 3D grid points, where Nq = x Ny x N^, Hq and 'Voi't) are Nq x Nq 
matrices with the matrix elements constructed by the unperturbed Hamiltonian ffo(k, k') and 
potential VD(k, k',t) on the 3D grid points together with the weights of quadrature, respectively. 

The P-space electron wavefunctions at any time can be calculated from Eq. ([3]) if initial electron 
wavefunctions are given. At finite temperature, the electrons populate on the single electron states 
of the unperturbed system according to the Fermi distribution function. While in the case of low 
temperature limit considered here, the electrons will occupy all the single electron states from the 
lowest energy up to the Fermi energy. Thus, the initial electron wavefunctions are the eigenfunctions 
of the unperturbed system below the Fermi energy. The eigenvalue equation for the eigenfunctions 
is obtained by setting t = in Eq. ([3]) and replacing the term on the left hand side of Eq. ([3]) by 
E^{0). Note that Vnit) = when t = 0. 

Applying the second-order split-operator method [25] to Eq. one obtains the propagation 
equation $(t + At) = Fd (r) PqPd (t) where. At is the time step size, t = t + At/2, Pq = 
g-iHoAt j^j^j Pd (r) = e-^Vi5W^*/2 are the propagators, and Vz)(t) = [Vd (t) + Vd (t + At)] /2. 
The propagators can be calculated using the eigendecomposition scheme of matrix [26]. If the 
uniform time step size (At = const.) is employed, the propagator Pq only needs to be calculated 
once. In contrast, the time-dependent propagator P/j (r) has to be calculated at every time step. 

Using the 7^-space electron wavefunctions, the current through the molecular device is calculated 

by 



/(t) = 3-^ J2 / dan. -Re / (ik(/)*(k, t)e 

(^^) j=occ ^c: 



-ikr 



dk'(l)j{k',t)k'e 



f Jk' r 



(4) 
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where, Sc is the surface enclosmg the C zone, n is the unit vector perpendicular to the surface 
element da, 0j(k, t) is the jth electron wavefunction with the momentum k, and the sum j is over 
all the occupied electron states of the C zone. 

To demonstrate the feasibility of the proposed P-space approach, we first apply it to a one- 
dimensional (ID) system driven by a DC bias. The potential of the unperturbed ID system is 
zero everywhere in the whole TZ space. The C zone extends from x = —6 to +6 a.u. Initially, all 
the single electron states are occupied up to the Fermi energy ep = 0.3 a.u. At t = 0, a DC bias 
Vh is applied to the electrodes. We choose vl = —vr = Vb/2, femax = 2.0 a.u., and At = 0.1 a.u. 



The identical ID problem has been considered in Ref. 



12, 



13] by using the transparent-boundary 



approach in the TZ space. In FIG. 1, we plot the current densities at x = as a function of time 
for three biases. It is shown that a steady-state current is achieved for each bias and approaches 
the steady-state current value obtained from the Landauer formula [l^ . [l^. Furthermore, the 
currents of the proposed 7^-space approach are in very excellent agreement with those of the TZ- 
space transparent-boundary approach [l^.fl^. However, the computational effort in the calculation 
with the T'-space approach is considerably smaller than that with the 7^-space approach. 

The second ID system we consider is a double square barrier potential driven by a time- 
dependent bias. The potential of the unperturbed system is v{x) = 0.5 a.u. for 5 < |x| < 6 
a.u. and zero otherwise [l^ . [l3 |. The C zone extends from x = —6 to -|-6 a.u. A time-dependent 
bias, Vb (t) = Vbo [6 {t) — 6 {t — to)], which is turned on at t = and turned off at t = 75 a.u., is 
applied to the electrodes, where 6 is the step function and to = 75 a.u. We choose ep = 0.3 a.u., 
Vl = Vb (t), vr = 0, /cmax = 2.0 a.u., and At = 0.1 a.u. In FIG. 2 we display the current densities 
at X = as a function of time for four biases. It is shown that for each bias the current oscillates 
after the bias is turned off and the amplitude is proportional to the bias. It is also shown that the 
currents tend to zero steady-state current quickly after turning off the bias and the transient time 
for turning off a bias is much shorter than that for turning on the bias [l^ . [l^ . Moreover, the 
results of the "P-spac e appro ach are in very good agreement with those of the 7^-space transparent- 



boundary approach 



12, 



131 ]. Again the proposed P-space approach reproduces the results of the 



7^-space transparent-boundary approach with considerably less computational effort. 

In summary, we propose a first-principles TDDFT approach in the V space for the study of elec- 
tron transport dynamics in molecular devices under the arbitrary biases. This approach is effective 
in theoretical treatment and efficient and feasible in computation. In this approach, the electron 
wavefunction is calculated by solving a time-dependent integrodifferential equation in a finite V- 
space volume. This equation is obtained by the Fourier transform of the 7?.-space TDKS equation. 
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It is exact and contains all the effects and information of the electron transport of molecular de- 
vices. It is free of the tricky self-energy function and memory term and beyond the extensively used 
WBL. In addition, unlike the 7^-space approach, the procedure based on this approach does not 
need to impose any complicated boundary condition to the electron wavefunctions. The computer 
resources such as CPU time used at each time step and the RAM required in the calculation with 
this approach do not increase with the number of time steps. The proposed P-space approach has 
been successfully applied to the study of ID systems with less computational effort, demonstrat- 
ing it promising to extend the proposed approach to the study of electron transport dynamics in 
realistic molecular devices 
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Biosciences Division of the Office of Basic Energy Sciences, Office of Science, U. S. Department of 
Energy, and by the National Science Foundation. 
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Figure captions 

FIG. 1 (Color online) Current densities at x = for the biases Vb = 0.1, 0.3, and 0.5 a.u. For 
each bias, the short dash curve is the result of the 7^-space transparent-boundary approach Q, [3], 
the dash line is the steady-state current from the Landauer formula [l^ , and the solid curve is 
the result of the 7^-space approach. 

FIG. 2 (Color online) Current densities at x = for the biases with v^q = 0.05, 0.15, 0.25 and 
0.35 a.u. For each bias, the short dash curve is the result of the 7^-space transparent-boundary 
approach [l^, [l^ and the solid curve is the result of the 7^-space approach. 
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